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ABSTRACT 



The problem of accurately replicating the parameters which 
define a given system for the purposes of implementing modern 
control strategies is important. Using an Autoregressive- 
Moving Average (ARMA) representation for the unknown system, 
a model is identified by processing input/output data to esti- 
mate the coefficients associated with the ARMA equation. Iden- 
tification of unknown system parameters using Kalman filtering 
methods was accomplished by augmenting the state vector. In 
this thesis the Kalman filter is formulated so that parameters 
can be identified explicitly. We call this approach the 
Adaptive Kalman Identifier (AKI) . 

It is shown that the Adaptive least mean square (LMS)f 
Adaptive Recursive LMS and Adaptive Lattice filters are special 
suboptimal cases of the AKI . The convergence and modeling 
properties are compared with those of the AKI by simulation 
using various types of data. 

With minor modifications, the AKI algorithm was used to 
identify the linear and non-linear ARMA models of the phase 
locked loop (PLL) . A discrete PLL using a forward Euler inte- 
gration scheme was used as a source of non-linear data. The 
AKI technique appears to enable one to discern when a potential 
non-linear system enters into its non-linear mode of operation. 
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I. 



INTRODUCTION 



The accomplishments in the area of microprocessor tech- 
nology in the last decade have made a noticable impact in the 
area of modern control. The desire to implement modern con- 
trol theories taking advantage of these advancements, has 
made it imperative that the engineer attempt to mathemat- 
ically replicate the system he ultimately desires to control. 
Efforts in this regard have, hence, generated a growing 
interest in the area of system identification [Ref. 1,2, 3, 4]. 

By implication this thesis concerns itself with discrete/ 
digital signal processing. 

A. BACKGROUND 

System identification or modeling can be accomplished by 
innovative application of existing techniques which were 
generally considered filtering or state estimation methods 
[Ref. 5], Previous research efforts have, for obvious reasons, 
focused on linear modeling; however, there is a rising interest 
in non-linear modeling methods [Ref. 6,7,8] . 

An attractive form for the representation of an unknown 
system is the Autoregressive -Moving Average (ARMA) equation, 

00 00 

y(k) = l a.u(k-i) - £ b.y(k-i) (1.1) 

i=0 1 i=l 1 

which states that the present output, y(k) , is a linear com- 
bination of past outputs, y(k-i) , and of past and present 
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inputs, u(k). Its attractiveness lies in its linear charac- 
ter and easily implementable structure using microprocessor/ 
computer algorithms. Its relationship to the state space 
representation of a linear system has already been established 
[Ref. 9,10] making it germane to consider that a system is 
identifiable by an equation of the form (1.1). The system 
identification problem thus entails identifying the coeffi- 
cients a . and b . . 

1 l 

There are several methods for computing the coefficients, 
a^ and b^, however, it is not the intent of this brief intro- 
duction to attempt to develop even the majority of them. 
Nevertheless, it is practical to present a referenced history 
of those methods encountered in this thesis. 

Adaptive algorithms for the purpose of estimating the 
coefficients of (1.1) have always been of interest. Widrow, 
using a least means square error criterion and implementing 
the method of steepest descent, developed an adaptive algorithm 
which estimated the coefficients of the moving average process 
associated with equation (1.1) [Ref. 11]. That is, the moving 
average model. 



y (k) 



a^u(k) + aju(k-l) +a^u(k-2) + ... 



( 1 . 2 ) 



where 




a 



0 
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= (a. 



' a O b 2 ) ‘ 



b l U l 



' a O b l> 



was identified. The LMS theory has since been extended to 
include block LMS filtering methods [Ref. 12,13,14] using 
various search techniques [Ref. 15,4(Chapter 5)]. These 
methods have enjoyed much popularity in the area of Linear 
Prediction and digital speech processing [Ref. 16,17]. 

The shortcoming of representing equation (1.1) by its 
equivalent moving average model (1.2) lies in the practical 
aspect of its implementation. That is, the infinite series 
represented by equation (1.2) must be truncated at some point 
resulting in an approximation which may not adequately repre- 
sent equation (1.1). Hence, efforts to adaptively estimate 
the coefficients for both the autoregressive (b^) and moving 
average (a^) processes continued [Ref. 18,19,20] with various 
degrees of success. Feintuch's Adaptive Regressive LMS pro- 
cedure [Ref. 18] is still a controversial issue [Ref. 19,20]. 

From yet a different direction, Anderson and Moore suggested 
that the Kalman filtering algorithms can be used as a means 
to identify the a^ and b^ coef f icients'*’ of equation (1.1) 

[Ref. 21, pp . 50-52]. This thesis exploits this application. 



Anderson and Moore in fact formulate the technique to 
compute the coefficients a. ,b. i = 1,2, ... p where 
p = n+m. 1 1 
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B. INVESTIGATIONS AND CONTRIBUTIONS 



This thesis formulates an adaptive application of the Kal- 
man filter to identify the coefficients of the general ARMA 
equation (1.1) . It is shown that the LMS adaptive filter 
[Ref. 11] and the Adaptive Recursive LMS filter [Ref. 18] 
are (1) special cases of the Adaptive ARMA Kalman identifier 
and (2) sub-optimal with respect to the underlying least 
means square (LMS) error criterion upon which the LMS adap- 
tive filter, the adaptive recursive filter and the Kalman 
filter are based. It is also demonstrated that the adaptive 
Kalman identifier has excellent convergence properties. 

By comparison, it is noted that the Kalman algorithm 
accounts for measurement noise where the LMS algorithms do 
not, a heretofore unapproached problem by LMS algorithms. 

The results indicate that the suggested modification by 
Griffiths [Ref. 22] of the convergence factor, k g , for the 
LMS adaptive filter is justified. 

The application of the Kalman filter algorithm is ex- 
tended to identification of the coefficients associated with 
a special case of the generalized non-linear ARMA model [Ref. 
8] . The results of the non-linear simulations suggest a tech- 
nique for determining when a potential non-linear system 
enters its non-linear operating regime . Such a technique can 
prove valuable when on-line performance evaluation of a known 
non-linear system is required. 

Lastly, the connection between the Kalman filter algorithm 
and the lattice filter algorithm [Ref. 6] is made. An example 



16 



which demonstrates and compares the performance of both 
algorithms is given. 

C. ORGANIZATION 

Chapter II presents and exploits the Kalman filter equa- 
tions emphasizing its connections with the Yule-Walker and 
the discrete Weiner-Hopf equations. The theory is further 
developed to investigate the general ARMA case. Chapter III 
pursues the theoretical comparisons between the discrete 
Weiner-Hopf equation upon which the adaptive LMS filter is 
based and the MA form of the Kalman filter. The theoretical 
comparison between the Adaptive Recursive LMS filter and the 
general ARMA form of the Kalman filter is made in Chapter IV. 
The software methods by which linear and non-linear synthetic 
data is generated are discussed in Chapter V. A short user's 
description of the data processing programs and the options 
provided is given in Chapter VI. A non-linear application 
for the identification of the parameters of a non-linear plant 
is presented in Chapter VII followed by a discussion and pre- 
sentation of the results in Chapter VIII. 
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II . 



FORMULATION OF THE PROBLEM 



Identifying the coefficients of the Autoregressive- 



Moving Average (ARMA) equation, 



y(k) + b 1 y(k-l) + ... + fc> n y (k-n) 



a^u(k) + a^u(k-l) 



+ ... + a u(k-m) (2.1a) 
m 



m n 



y (k) = l a.u(k-i) - £ b.y(k-i) 

i=0 1 i=l 1 



(2.1b) 



can be formulated as an adaptive Kalman identification problem. 
That is, instead of using the well known Kalman Filter equa- 
tions [Ref. 23,24] to recursively estimate the states of a 
system, one can utilize the Kalman filter equations to adap- 
tively estimate the coefficients of either an Autoregressive 
(AR) , Moving Average (MA) , or Autoregress ive -Moving Average 
(ARMA) process by proper definition of the quantities involved. 
We call this the Adaptive Kalman Identifier for obvious 
reasons . 

A. THE DETERMINISTIC CASE 

Consider for the moment that the a^ and b^ are constant. 
Then by collecting a sufficient number of measurements of the 
input u(k) and the output y(k) , one can readily solve a set 
of linear equations for the a^ and b,. coefficients. The 
"sufficient number" that is needed is n+m+1 which define the 
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n+m+1 linear equations which solve the n+m+1 unknown coeffi- 
cients. The matrix equation to be solved takes the form: 



y(0) 




u(0) 


0 


... o ! o 

i 


• • • 





0 




a o 


yd) 

• 




u(l) 


u(0) 


! y(0) 

1 

1 

1 






0 




a 

1 


• 

y(m) 


— 


u(m) 


u (m-1) 


• t ft • 

• • • 

i 

u (0) J y (m-1) 

i ■ - 


y (m-2) 


y(0) 


• 

• 

0 




a m 


y(m+l) 




u(m+l) 




u(l) j y (m) 

I 


y (m-1) 


y(l) 


0 




b i 


• 








1 

1 

1 

1 

1 

1 










b 2 


y (m+n) 




u(m+n) 


• • • 


1 

I 

u(m) j y(n-l) 


y(n-2) 




y(0) 




-1 



( 2 . 2 ) 



Z 




(2.3) 



It is interesting to note that the (n+m+1) x (n+m+1) matrix, 
T, is block symmetric. The solution to (2.3) is readily 
apparent, namely. 



= T -1 y; (2.4) 

where T ^ is the inverse matrix of T. Since we have assumed 
that the elements of T are perfect noiseless measurements, 
and that we somehow know a-priori the number of unknown 



a 

b 
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coefficients, then T is of full rank, namely rank (T) 
n+m+1, and the inverse of T exists. 



B. THE NONDETERMINISTIC CASE 

Rarely can the coefficients be modeled so ideally. A 
more prudent and realistic model admits that the ARMA equa- 
tion coefficients are subject to random perturbations. Further, 
it can be said that in general measurement devices introduce 
noise into the measurement data. Hence, the measurement 
model should reflect this fact. Developing the Adaptive Kalman 
Identifier along the same lines as [Ref. 21] we let. 



a^ (k+1) = a^ (k) + w^ (k) i = 0,1,2, ... m (2.5a) 

bj (k+1) = b_.(k) + Wj (k) j = 1,2, ... n (2.5b) 

where the (w^ (k) ,w^ (k) } are samples from a zero mean, white, 
gaussian random process. Additionally, we assume that the 
noise sources are uncorrelated, 



E [w (k) w (k) ] 
a r a s 


= 0 


for 


r ^ s 


(2.6a) 


E (w (k)w (k) ] 
r s 


= 0 


for 


r ? s 


(2.6b) 


E [w (k) w (k) ] 
a r b s 


= 0 


for 


all r,s 


(2.6c) 


Similarly, allowing 


for noisy measurement devices , 


the 


measurement equation 


(2.1b) 


, is 


modelled as, 




-H 

1 

3 

♦H 

fd 

o 

II 

•H 

II 

>i 


n 

• l 

i=l 


b^y(k-i) + v(k) 


(2.7) 
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where the {v(k)} are samples from a zero mean, white, gaussian 
random process. It is also assumed that, 



E[w i (k) v (k) ] 



0 for all i 



(2.8a) 



E[ Wj (k) v ( k ) ] 



0 for all j 



(2.8b) 



That is, it is assumed that the noise perturbations associated 
with one coefficient are independent of the noise perturba- 
tions associated with any other coefficient, and that the 
measurement noise and the coefficient perturbations are 
independent. 

At this point it is ncessary to use judgment and experi- 
ence and utilize all the information known about the physical 
system which equation (2.7) represents to assign variances for 
the random processes (w^(k)}, (w^ (k) } and (v(k)}. Let 



vsr (k) 



Q = E 



[w i (k) Wj (k) ] > 



( 2 . 9a) 





0 



Q 



(2.9b) 



0 




2 

where, Q 1 = diag ( a ) , Q 2 



i 



diag (a^ ) and, 
r w . 

J 



R = E [v (k) v (k ) ] 



(2.9c) 



and E[x] is the expected or average value of x. 
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Instead of using x's denoting the commonly used notation 
for state variables, we retain the flavor of the problem by 
using the a^'s and b^ ' s as the "states" of the Adaptive Kalman 
Identifier. It is hoped that a more meaningful understanding 
of the Adaptive Kalman Identifier may thus be gained. There- 
fore, define the state vector. 



a 0 (k) 

a x (k) 



a (k) 
m 




b 1 (k) 

b 2 (k) 




( 2 . 10 ) 



b n (k) 

Combining equations (2.5) and (2.10), we have the discrete, 
first order Gauss-Markov process [Ref. 3], 



a (k+1) 



b (k+1) 



a (k) 
b(k) 



+ w (k) 



whe re , 



w (k) 



w i (k) 



w . (k) 
“3 



( 2 . 11 ) 



( 2 . 12 ) 
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Note that w(k) is an n+m+1 vector whose elements are a con- 
catenation of the noise sequences associated with the a i and 
b^. To complete the Adaptive Kalman Identifier formulation 
we define the measurement vector, 



H(k) = [u(k) u(k-l) ... u(k-m) -y(k-l) ... -y(k-n)] (2.13) 



and its associated measurement equation, 



y (k) = H (k ) 



a (k) 
b (k) 



+ v ( k ) 



(2.14) 



Then the solution to the Adaptive Kalman Identifier problem 
[Ref. 21] is, 



a (k+1 Ik) 



b (k+1 |k) 



= [I -K (k)H (k) ] 



a(k Ik-1) 



b (k I k-1) 



+ K(k)y(k) (2.15a) 



K (k) = P (k I k ( 1 ) H T (k) [H (k) P(klk-l) H T (k) +R]" 1 



(2.15b) 



P (k+1 1 k ) = P (k | k— 1 ) - K(k)H(k) P(k|k-1) + Q. 



(2.15c) 



Equations (2.15) are initialized by assuming an initial value 
for the coefficients (2.10) and assigning to our assumption 
a measure of our confidence in the initial guess. That is, 
we pick the values, 



a ( 0 | -1 ) 



b (0 | -1) 



and P (0 -1) 



(2.16) 
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where P (0 -1) is defined as the error covariance of the 
coefficients , 



P (k k-1) 




a (k |k) 
b (k | k) 




a(k) 
b (k) 



(2.17) 



for k = 0. Equation (2.15b) is generally referred to as the 
Kalman gain. Two special cases are of interest: case (1); 
a. = 0 for all i, and case (2); b. = 0 for all i. 

l l 

1. Autoregressive Form (a^ = 0) 



For case 1, equation (2.7) takes the form, 



n 

y(k) = - l b.y(k-i) + v(k). (2.18) 

i=l 1 

This is a recursive equation stating that the present output 
is a linear combination of past outputs corrupted with addi- 
tive gaussian noise, v(k). Equation (2.18) is more formally 
recognized as an autoregressive prcoess of order n [Ref. 25, 
26] . Writing the Kalman solution for the coefficients of the 
AR process. 



b(k+l|k) = [I - K (k) K (k) ] b (k | k-1 ) +K(k)y(k) (2.19a) 

K (k) = P (k | k — 1 ) H T (k ) [H (k)P(k |k-l) H T (k) + R] -1 (2.19b) 
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P (k+1 | k) = P (k | k-1) - K(k)H(k)P(k|k-l) + Q. (2.19c) 

Alternate forms of the Kalman equations given by Maybeck 
[Ref. 27] are, 

b (k+1 1 k) = P (k | k-1) P (k | k-1) -1 b(k|k-l) 

+ P (k+1 |k) H T (k) R -1 y (k) (2.20a) 

P (k+1 |k) = (P (k | k-1 ) ) _1 +H T (k) R -1 H(k) _1 . (2.20b) 

We can model the fact that we have no a-priori knowledge 
about the initial values of the coefficients by letting, 

[P (0 [ -1 ) ] _1 ~ 0. (2.21) 

Further, if we remain ignorant and totally doubt our previous 
estimate, then P(k|k-1) is modeled as, 

[P (k | k— 1 ) ] -1 ^ 0. (2.22) 

Equations (2.20) reduce to, 

b (k+1 |k) = [H (k ) R -1 H T (k)] _1 H(k) R _1 Y(k), (2.23) 

the weighted least squares estimate, previously encountered 
by Swerling [Ref. 28,29,30], of the coefficients, b. Carry- 
ing the analysis further and letting R ^ = 1, we have the 
Penrose pseudo-inverse solution [Ref. 31,32], 

b (k+1 | k) = [H (k) H T (k) ] -1 H (k) y (k) (2.24) 
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Swerling [Ref. 28] has shown that if the weighting matrix 

in (2.23) is not the inverse of the covariance matrix of 

the measurement errors , then the accuracy of the estimated 

coefficients (2.24) will be degraded. 

The aforementioned notwithstanding, we press further 

T 

into the analysis of equation (2.24) . The product H(k)H (k) 
can be written as, 



y (k-1) 
y (k-2) 



H (k) E (k) = 



y (k-n) 



[ y (k-1) y (k-2 ) ... y(k-n)] 



H(k)H (k) = 



y (k-1) y (k-l)y (k-2) ... y(k-l)y(k-n) 

y (k-2) y (k-1) y 2 (k-2) 



y (k-n) y (k-1) 



. y (k-n) y (k-n) 



(2.25) 



(2.26) 



As we let k ■* <*> the expected value of (2.26) becomes the 
autocorrelation matrix. 



E{ lim H (k) H T (k) } = R (k) 

k+- yy 



(2.27) 



where , 
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(2.28) 





Similarly, the product H(k)y(k) as k -*■ °° becomes 





(2.29) 




The steady state solution for the estimate of the coeffi- 



Equation (2.30) is one of the starting points from which 
Perry [Ref. 6] develops his Lattice modeling algorithms. 
Equation (2.30) can also be recognized as the solution to 
the Yule-Walker or Normal equations [Ref. 26]. 

2. Moving Average Form (b^ = 0) 

Referring once again to equation (2.7) and setting 
the coefficients b^ = 0 for all i we have case (2) , 



cients b(k+l|k) is. 



b (k+1 Ik) = R -1 (k) r (k) . 

— ss 1 yy yy 



(2.30) 



m 



y (k) = l u(k-i) + v (k) 

i=0 



(2.31) 
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Since v(k) is by definition noise associated with the 
measurement, we can combine y(k) and v(k) such that, 

m 

z(k) = y (k) - v (k) = £ a.u(k-i). (2.32) 

i=0 1 

Equation (2.32) simply states that the present measurement 
is a linear combination of past and present inputs or by 
definition, a Moving Average (MA) process [Ref. 25,33]. The 
Adaptive Kalman Identifier estimate for the coefficients of 
the MA process is as before, 



a (k+1 |k) = [I - K (k) H (k) ] a(k+l|k) +K(k)z(k) (2.33a) 

K (k) = P (k | k-1) H T (k) [H (k) P (k | k-1) H T (k)+R] -1 (2.33b) 

P (k+1 |k) = P (k ] k-1) - K(k)H (k)P (k | k-1 ) +Q. (2.33c) 

Using the alternate forms of the above equations we arrive 
at. 



a (k+1 1 k) = [P (k+1 1 k ) [ P (k |k-l ) ] x ]a(k|k-l) 

+ [P (k+1 | k) H T (k) R -1 ] z (k) (2.34a) 

P (k+1 |k) = [ [P (k 1 k-1 ) ] _1 + H T (k) R _1 H(k)]" 1 . (2.34b) 

Arguing as we did for the autoregressive case, no a-priori 
knowledge about the initial values of the moving average 
coefficients, implies that, 
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0 . 



(2.35) 



[P(0 -1) ] 



-1 



And even if after we processed one measurement we still 
admitted no knowledge as to the accuracy of the previous 
estimate, we imply that. 



[P(kik-l)] 



-1 



(2.36) 



Substituting these implications into equations (2.34) our 
estimate becomes. 



a (k+1 Ik) = [H (k) R _1 H T (k)] _1 H(k) R _1 z(k) 



(2.37) 



the by now familiar weighted least squares estimate [Ref. 28, 
29,30]. If once again we allow R ^ to be unity, let k tend 
toward infinity and take the expectation of equation (2.37), 
one arrives at the discrete form of the Wiener-Hopf equation 
[Ref. 11] , namely 



a (k+1 Ik) = R -1 (k) r (k) 
— 1 uu uz 



(2.38) 



where it can easily be shown that. 



R (k) = 

uu 



R (0) 
uu 



R (~m) 
uu 



R (-1) 
uu 



R (-1) R (0) 
uu uu 



R (-m) 
uu 

R ( 1-m) 
uu 



R (0) 
uu 



(2.39) 



and, 
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